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Abstract. Lattice Boltzmann methods are numerical schemes derived as a kinetic approximation of an underlying lattice 
gas. A numerical convergence theory for nonlinear convective-diffusive lattice Boltzmann methods is established. Convergence, 
consistency, and stability are defined through truncated Hilbert expansions. In this setting it is shown that consistency and 
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are presented. 
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1 Introduction 

Lattice gases, which were introduced in the early 1970s [14,15], have been used to simulate problems in 
fluid dynamics [8, 11, 12]. A Lattice gas involves indistinguishable pseudo-particles that traverse from node 
to node along the links of a lattice in unison according to the ticks of a discrete clock and that interact 
at the nodes of the lattice. An exclusion principle is imposed so that the state at any given node may 
be described with a finite number of bits. Thus, lattice gases are amenable to a mathematical description 
over a Boolean field and have been related to cellular automata [22,25]. The microscopic evolution of a 
lattice gas system can be viewed as a space-time- velocity discretization of the Boltzmann equation (see, e.g., 
[2]), in which the precision of the particle distributions is reduced to as few as one bit. Characteristic of 
lattice gas methods is that the velocity discretization remains fixed (to the lattice structure) in the limit 
as the spatial and temporal discretization parameters tend toward zero. While the microdynamics of a 
lattice gas is certainly not physical, the aim is however to recover physical macrodynamics via this simple, 
non-physical microdynamic means. This and the accompanying statistics have been explored for a variety 
of lattice gas methods [8] . In lattice Boltzmann methods particle distributions (not particles) traverse the 
links of the lattice and interact at the nodes thereof [18, 20, 24] (see, e.g., [8], for further references). Lattice 
Boltzmann methods do not possess the statistical fluctuations that are inherent in lattice gas methods. For 
certain classes of these methods we develop a convergence theory. Such classes include linear and nonlinear 
convective-diffusive and monotone lattice Boltzmann methods. 

Our paper is organized as follows. In the next section we quantify the microscopic dynamics of lattice 
gases and then derive formally an equation for the expected or mean behavior of this system, the so-called 
lattice Boltzmann equation, which forms the basis of lattice Boltzmann methods. This equation has the form 
of a discrete space-time kinetic equation composed of an advection part and a collision part, the so-called 
Boltzmann collision operator, whose properties are examined in section 3. The main thrust of this paper is to 
establish a convergence theory for solutions of this equation in the continuum limit for a class of convective- 
diffusive lattice Boltzmann methods. In section 4, we identify the class of lattice Boltzmann methods that 
have a convective-diffusive continuum limit through an analogue of the classical Hilbert expansion of kinetic 
theory. This is the lattice Boltzmann equivalent of the consistency step of traditional convergence proofs 
for numerical schemes. Stability, and therefore convergence, is then established in section 5 for a class of 
so-called monotone lattice Boltzmann methods. Specific examples of both diffusive and convective-diffusive 
lattice Boltzmann methods that are both consistent and monotone are presented in sections 6 and 7. 

2 Lattice Gas Dynamics 

A lattice gas [8] involves indistinguishable particles moving about from node to node on a lattice in unison 
with the ticks of a discrete clock and interacting at the nodes of the lattice. More precisely, each particle is 
characterized as being in one of a finite number of possible particle states, and associated with each possible 
state is a velocity, which is the lattice vector on which the particle will translate during the advection step 
of each clock cycle. Before the advection step of each cycle, the gas undergoes a collision step during which 
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the particles at each node interact according to a set of rules that change the states of the particles at that 
node independently of the state of the particles at any of the other nodes. These collision rules may be either 
deterministic or stochastic. 

A lattice gas automaton [25] envokes the exclusion principle, which states that at any given time there is 
at most one particle per particle state per node. This principle ensures that a single bit, called an occupation 
number, can encode the absence (=0) or presence (=1) of a particle in a particular particle state at a given 
node. Thus, a finite number of bits may be used to describe the local state of the gas at a given node — one 
bit for each possible particle state. 

More specifically, a lattice gas automaton and its dynamics are quantified as follows [11]: 

1. A spatial lattice domain, X C H D . More precisely, given a macroscopic domain ft C R D , set 
X = f2 fl L for some regular Z?-dimensional lattice L C R D with microscopic spacing Ax. The nodes 
of X are denoted i. In order to avoid the complications wrought by boundaries, we shall assume that 
X and O are effectively without boundary by imposing periodicity, say of length L. 

2. The ticks of the discrete clock are called cycles and are indexed by to = 0,1,2,.... Each cycle 
corresponds to a microscopic time step of At. 

3. A finite set P (of cardinality |P|) of possible particle states at each node. A mapping v : P — > L 
associates a velocity vector v(p) with values in a lattice neighborhood of the origin to each particle 
state peP. In the absence of collisions, a particle in state p will translate each cycle by v(p) along 
the lattice. 

4. The absence (=0) or presence (=1) of a particle in a particular particle state p at the node i after 
cycle m is encoded by an occupation number N(m,i,p) G 2 = {0, 1}. The local state at the node 
i after cycle to is given by N(m, i, •) G 2 P . (A "•" in an argument denotes a function over the dotted 
argument.) 

5. The advection operator A translates the particles to neighboring nodes and advances the discrete 
time cycle from to to to + 1. It is defined by 

AN(m,i,p) = N(m+l,i + v(p),p) . 

6. The collision operator C acts locally in space-time lattice to determine the change in local state due 
to interactions between particles. More specifically, given the local state N(m, i, •) at node i after cycle 
to, the collision operator determines a collided state TV' (to, i, •) by 

AT' (to, i, •) = N(m, i, •) + C{N(m, i, •)) • 

In other words, the map n' = n + C(n) takes the set 2 P of local states into itself. 

7. The composition of the advection step with the collision step gives the microdynamical equation 

for cycle to + 1 as 

AN(m, i,p) = N(m, i,p) +C(N(m, i, -))(p) , 
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or, more simply, 

(1) AN = N + C(N). 

This equation states that after cycle m + 1 the new occupation number for state p at the new location 
i + v(p) is the same as the occupation number for state p at the location i after cycle m, plus some 
collisional correction. 

Since N(m,i, •) e 2 P , it can take on one of 2' p l possible local states in 2 P . In order to detect exactly 
which state is occupied, the general expression for the collision operator requires a representation of the 
Kronecker delta. Let N, n e 2 P have components N(p),n(p) and define 

N n — Y[ N(p) n ^ , N n = Y[ (1 A r (p)) (1_ " (p)) • 

Notice that N n N n = d(N,n), the Kronecker delta. The collision rules of the lattice gas determine a unique 
post-collisional state n' for any given prccollisional state n e 2 P . Introduce a matrix S(n 7 n') such that 



S{n,n') = 



1 , 

, otherwise . 



The corresponding collision operator can then be expressed as 

C(N)= (n -n)N n N n S{n,n'). 

n,n'e2 p 

In general, the matrix S may depend on node and/or cycle, even when the gas is deterministic. For example, 
it may take on alternate values at odd and even cycles, or at adjacent nodes, or both. Since S represents 
the collision map (every state n must have exactly one image n'), it satisfies 

(2) n ') = 1 . for ever y ne2 p . 

n'£2 p 

It is also clear that the collision map is one-to-one if and only if S satisfies 

(3) <^( n ' n ') = 1 ' for ever y n ' e 2P 

For simplicity, we will model all lattice gases as time stationary, spatially homogeneous stochastic pro- 
cesses. Let <S be the expected value of S. Then 



S(n, n) 



1 , with probability S(n, n!) ; 
, with probability 1 — S(n, n') . 



Of course, if the gas has only one possible collision map S, then S = S. The 2l p l x 2l p l matrix S is called 
the local transition matrix of the lattice gas. By (2), S satisfies 

(4) ^ S{n,n') = 1, for every n e 2 P . 
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It is also clear from (3) that if every possible collision map is one-to-one then S also satisfies 

(5) n ') = 1 ' for cvcr y n ' e 2P - 

n£2 p 

Let F — ExpVal(iV) where N is determined from the microdynamical equation (1). Then F takes on its 
value in the |P|-dimensional interval 1 = [0, l] p . Wc consider the equation 

(6) AF = F + C(F), 
where 

(7) C(F)= (n 1 -n)F n F rn S{n,n'). 

n,n'e2 p 

Here C is called a Boltzmann collision operator. If the expected value operation passes through the 
nonlinearities of the collision operator C, of the lattice gas, then F = F. Clearly, A(F) = F + C(F) g I. 
Equation (6) can be viewed as a finite difference equation whose solutions are grid functions F(m, i, •) where 
F(0, i, •) = i*b(i, •) is the initial condition. This equation is called the lattice Boltzmann equation for the 
lattice gas automaton [4-7,9-17]. 



3 Boltzmann Collision Operators 

We examine the relationship between the concepts of conservation, equilibria and dissipation for a Boltzmann 
collision operator C given by (7). These relations are not special to convective-diffusive lattice gases [8], but 
rather very general. The discussion here emphasizes this generality and closely parallels the treatment in [1] 
of Boltzmann collision operators without exclusion terms. 

The sum of any scalar or vector valued function / = f(p) over the variable p will be denoted by (/): 

</> = £/(?)■ 

p€P 

Concepts of conservation arc central to the existence of macroscopic limits. Two of them appear below 
which will be shown to be equivalent for collision operators satisfying suitable conditions. 

Definition 3.1 A mapping e : P — > R (alternatively a vector e g R p ) is said to be a locally conserved 

quantity for the collision operator C if (eC(F)) = 0, for every F el. 

Note that a locally conserved quantity e for the operator C given by (7) satisfies 

J2 {{n 1 -n)e)F n F rn S{n,n') = Q, 

n,n'e2 p 

for every F g I. Since the family of polynomials parameterized by n and given by f h p n p n [ s linearly 
independent, the above equality holds if and only if the coefficient of each F n F n vanishes: 

(8) J2 W -n)e)S(n,n')=0, 

n'e2 p 

for every n g 2 P . 

Another notion of conservation is one that holds for individual collisions. 
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Definition 3.2 A vector e e R p is said to be a microscopically conserved quantity for the collision 
operator C given by ( 7) if 

(9) ((n' - n)e)S(n,n') = 0, 
for every n, n' £ 2 P . 

While it is clear from comparing (8) and (9) that any microscopically conserved quantity is also locally 
conserved, the converse is generally not true. 

The converse is true for the following class of collision operators, however. 

Definition 3.3 The operator C is said to be in detailed balance if 

(10) S(n' , n) = S(n,n r ) , for every n,n' S 2 P , 
and in semidetailed balance if (see (4) and (5)) 

(11) ^ S(n',n)= ^ S(n,n') = l, for every n e 2 P . 

n'e2 p n'e2 p 

Clearly, the notion of semidetailed balance is a weakening of the detailed balance condition. Its usefulness 
arises through the following characterization, the proof of which is immediate. 

Lemma 3.4 The collision operator C given by (7) is in semidetailed balance if and only if 

(12) Yl "(n)S(n,n')= ]T v(ri)S(n,ri), 

n,n'G2 p n,n'e2 p 

for every v : 2 P — > R. 

The first implication of semidetailed balance is the following. 

Theorem 3.5 For any collision operator C given by (7) that is in semidetailed balance, any locally conserved 
quantity is also microscopically conserved. 

Proof: Let e e R p be a locally conserved quantity. Multiplying (8) by (ne) and summing over n gives 

(13) = ^ ((n' -n)e)(ne)S(n,n') 

n,n'£2 p 

= — ((ne) 2 — (n'e)(ne)) S(n,n) . 

n,n'£2 p 

Applying (12) of Lemma 3.4 to half of the first term inside of the last sum of (13) gives 

= (|W 2 + |W 2 - (n'e)(ne))S(n,n') 

n,n'e2 p 

= Yl \((n-n')e) 2 S(n,n'). 

n,n'e2 p 
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Since each term of this last sum is nonnegative then all of them must be equal to zero. But that implies (9) 
is satisfied and shows that e is also microscopically conserved. □ 

The set of all locally conserved quantities of C is a linear subspace of R p denoted by E and assumed 
to be nontrivial. Let K be the dimension of E and {ej \j — 1, . . . , K} be a basis. Let e : P — > where 
the components of e{p) are these basis vectors. Vectors in R K will be denoted with arrows. The Euclidean 
inner product of two such vectors, a and /3, is denoted by a © (3. 

Boltzmann collision operators in scmidetailed balance also satisfy the following iJ-theorem. 

Theorem 3.6 (_ff-Theorem) Suppose the collision operator C given by (7) is in semidetailed balance. Then 
it has the dissipation property 

(14) (C(F)Iog(F/F))<0, 

for every F £ I. Moreover, the following characterizations of equilibrium are equivalent: 

(i) C(F) = 0, 

(15) (ii) (C(F)log(F/F))=0, 

(Hi) F n F n =F n Vn,n' £ 2 P such that S(n,n') > 0, 
(iv) F = M0) for some (3 £ R K , 

where M([3) is given by, 

(16) M0) = - ^— -. 

1 + exp(-/3 e) 

Here e is a basis o/E, the K -dimensional space of locally conserved quantities. 

Remark: The form of the local equilibrium given in (16) is that of a Fermi-Dirac density, which is the quan- 
tum mechanical analogue of the classical Maxwellian density for particles satisfying an exclusion principle, 
hence the designation M. 

Proof: Since the logarithm of a product is the sum of its logarithms, one can verify that 

pn' pn \ 



(17) ((n-n')log(F/F)) = -log 



pn pn j 

Since C is in semidetailed balance, letting v(n) = F n F n in (12) of Lemma 3.4 yields 

0= Yl F n ' F n ' S{n,n') - F n F n S(n,n') . 

n,n'e2 p n,n'e2 p 

Hence, 

(18) -(C(F)log(F/F)> = J2 ((n-n')log(F/F))F n F n S(n,n') 

n,n'£2 p 

= y - — £- 1 - log - — F n F n S{n, n') . 



n,n'e1 p 
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Since y — 1 — logy > for every y € R+, every term in the last sum of (18) is nonnegative, so that the 
collision operator C satisfies the dissipation property. 

The characterization of equilibria (15) is argued as follows: (i) implies (ii) implies (iii) implies (iv) implies 
(i). The first implication is obvious. Assuming (ii), the last sum in (18) is zero and each of its nonnegative 
terms must vanish. This gives the formula 

^pn'pn' _ pnpn _ pnpn bg ^ p n F n S{n, „')=(), 

for every n, n' G 2 P . Since (y, z) i— > y — z — z \og(y/z) over R + x R + is nonnegative, and vanishes only on 
the diagonal y = z, it then follows that 

(F n F n - F n 'F n ')S(n,n') = 0, 

for every n, n' e 2 P , which gives (iii). Assuming (iii) then (17) implies 

(pn' pn' \ 
-_-j5(n,n')=0, 

for every n,n' G 2 P . Thus log(F/F) satisfies (9) and is therefore a microscopically conserved quantity. 
Hence, 

K 

log(F/F) = J2P j e j =0@S 

3 = 1 

for some vector j3 = (/3i, . . . , (3k) T € R^. Solving this for F yields (iv). Finally, assuming (iv) and using the 
fact that all locally conserved quantities are microscopically conserved (Theorem 3.5) and employing (12) of 
Lemma 3.4, it is easy to show that C(M((3)) = for every (3 e K K . □ 

Finally, another consequence of the property of semidetailed balance is the characterization of the Fred- 
holm alternative for the first derivative of C evaluated at any given local equilibrium M — M({3). The 
linearized collision operator at M is C M defined by 

(19) C M h = VC(M)h = d e C(M + eh) 

for every h e X. First observe that every e € E can be written as e = a e for some unique a € R K . The 
formula for the local equilibria (16) then yields 

(20) a®dgM0) = M0)(1 - M0))a Q e = M0)(1 - M0))e . 
Defining W M = W M 0) = M0){1 - M0)), one sees that 

(21) = a® d C(M0)) = VC(M0))(ctQ d M0)) =VC(M0))W M 0)e , 

which implies that E C Null(£ M VF M ). 

Next, let £^ denote the transpose of C M defined by 

(22) (hClg) = (g£ M h) , for every h,g£l. 
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If e is a locally conserved quantity of C then 



(eC M h) = d e (eC(M + eh)} 



= 0, 

e=0 



for every del, which by (22) implies that E C Null(££). 

The above inclusions become equalities for collision operators that are in semidetailed balance. 

Theorem 3.7 If the collision operator C given by (7) is in semidetailed balance and M is any local equilib- 
rium of C, then its linearization C M defined by (19) satisfies 

E = Null(£ M W M ) = Null(^) = Null(£ M W M + W M C T M ) , 

and moreover, (gC M W M g) < for every g £ E. 

Proof: Above it was shown that E is contained in both Null(£ M W^f) and Null(£^). Since every g e R p 
satisfies 

(gC M W M g) = \{g{C M W M + W M C T M )g) = {gW M C T M g) , 

it is clear that Null(£^) and Null(£ M W M ) are each contained in Null(£ M W M + W M C"^). All that remains to 
be shown is that Nul\(C M W M + W M C^) c E. A direct calculation following (19) and using C(M) = yields 

(23) £ M W M g= Yl ( n ' -n)(ng)M n M n S(n,n'). 

n,n'£2 F 

If g € Null(£ M H^ + W M Cj 1 ) then using semidetailed balance (as in the proof of Theorem 3.5) shows that 
0=-(g£ M W M g) = - £ ((n' - n)g){n g)M n M n S(n, n') 

n,n'e2 p 

J2 \{{n' -n)g) 2 M n M rn S{n,n'). 

n,n'e2 F 

Since each term of this last sum is nonnegative, all of them must be equal to zero. But that means g satisfies 
(9) and is therefore a locally conserved quantity (g € E). Moreover, it is clear that the sum is zero if and 
only if g G E. □ 

An immediate consequence of Theorem 3.7 is the Fredholm alterative that for any / E R p the overdetermined 
system 

(24) £,/! = /, (eh)=0, 

has a solution if and only if (ef) = 0, in which case the solution is unique and is denoted by f . The 
operator L~} is a (left) pseudoinverse of C M . 



4 The Hilbert Expansion and Diffusion 

Here we give the characterization of convective-diffusive lattice gases by properties of their Boltzmann 
approximations, more precisely, by properties of their local equilibria. The notion of the continuum limit of 
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such a gas involves refining the lattice domain X within the macroscopic domain il E R D and is formulated 
in terms of the vanishing of a parameter S > that is related to lattice spacing Ax and time cycle interval 
At by 

(25) Ax = 5L, At = S 2 T, 

where L, T > are macroscopic length and time scales. Of course, the scaling of At = 0((Ax) 2 ) is the usual 
diffusive scaling, but not every lattice gas has macroscopic dynamics that is consistent with it. 
Here we consider Boltzmann collision operators of the form 

(26) C = C (0) +SC W , 

where and arc Boltzmann operators such that every locally conserved quantity of is also locally 
conserved by C^h The spaces of locally conserved quantities of and C therefore coincide and we denote 
this space by E and let e denote a basis. Moreover, we assume that C(°) is in scmidctailcd balance; hence its 
equilibria are given by F — M(/3), and it satisfies the if-Theorem. Finally, we assume that the lattice gas 
is diffusive: 

Definition 4.1 A lattice gas such as given above is called diffusive provided 

(27) (veM(/3)) = , for every [3 e H K . 

This condition will insure that the time scale of the macroscopic dynamics will be consistent with the diffusion 
scaling (25). 

The limiting convective-diffusive macroscopic dynamics of the gas is established as follows. First, a family 
of approximate solutions of the lattice Boltzmann equation parametrized by S is constructed from smooth 
functions over the (t, x)-domain R + x f2 that are solutions of convection-diffusion equations. Then it is 
shown that the exact solution of the lattice Boltzmann equation and the approximate solution converge in 
some sense. The first step, carried out below, is the lattice Boltzmann version of the consistency step of 
most numerical convergence proofs, while the second will follow from a stability argument given in the next 
section. 

Given any function H = H(t,x,p) that is a smooth mapping from the (t, x)-domain R + x f2 into R^, 
the Taylor expansion of AH(t,x,p) about (t, x) is 

AH(t,x,p) = H(t + 5 2 T,x + 5Lv(p),p) 

(28) = Y,^[S 2 Td t + SL V (p)-WyH(t,x,p). 

3=0 J ' 

Grouping terms by order of 5 gives 

AH-H= 5L{wV)H 

+ 5 2 [Td t + \L 2 (wV) 2 ]H 

(29) + S 3 [TLdt(vV) + ±L 3 (vVf]H 

+ 5 A [\T 2 d 2 + \TL 2 d t (w ■ V) 2 + ±L 4 (v ■ V) 4 ]H 
+ ■■■ . 
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We construct an approximate solution to the lattice Boltzmann equation 

(30) AH-H = C(H) 
by formally expanding H in powers of S as 

oo 

(31) iT(t,x,p) = £W fc >(t,x,p). 

k=0 

This series is the lattice analogue of the classical Hilbert expansion of kinetic theory through which one 
formally passes to the limiting macroscopic dynamics. Notice that, as with classical Hilbert expansions, the 
advection side (29) is O(S). 

Expanding the left side of (30) in powers of S gives 



(32) AH-H = (A- I)J2 skh(k) 



k=0 



where 



= a^+Sa^+S 2 a^+S 3 a^ + --- , 

a(°> = 0, 

a« - L(vV)/i< >, 

= [Td t + ±L 2 (vV) 2 ]h^ +L(vV)h {1 \ 
a< 3 ) = [TLd t {w • V) + ±L 3 (v • V) 3 ] 

+ [Td t + \L 2 (y • V) 2 ] + L(v • V)h™ , 



Notice that depends on through 1 \ but not on 
Similarly expanding the right side of (30) gives 

(33) C(F) = cw(f^S k h (k A+5C^(jr6 k h( k A 

^fc=0 ' ^fc=0 ' 

= c<°> + 5cW +5 2 cW +^ c (3) + ... 5 



where 



.(0) = c<°>(/i<°>), 

:« = VC(°\hW)-hW+CW(hW), 

= VC<S>\hW) ■ + \V 2 C^\h^) ■ hWhW+VCW{h<V) ■ , 

:( 3 ) = VC ia \h^)-h^ +V 2 C^\h^)-h^h^ 

+ VCW{hM) ■ + \V 2 C^\h^) ■ hWhW , 
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The general form for the in (33) is 

(34) =T>C<®(hM)-hW+rW , 

where refers to all the remaining terms of , each of which depends on through h^ -1 ^ , but not on 
h^K Notice that since is just the sum of derivatives of the collision operator C, it automatically satisfies 
(e>W) = 0. 

Matching (32) to (33) at leading order gives 

(35) C<®(h,M) = 0, 
which by the i?-Theorem implies that 

(36) h (0) = M = M0) , for some (3 E K K . 

Here (3 = (3(t, x) is a smooth function still to be determined. 

Matching (32) to (33) order by order for j > gives a linear equation for in the form 

(37) C M h^=a^-r^\ 

where the right side depends on through but not on h^K Being of the form (24), the linear 

equation (37) will have a solution if and only if its right side satisfies the solvability condition 

= (e>W-r( j) )) . 

Since automatically satisfies (er^) = 0, this condition reduces to 

(38) = . 
This satisfied, the general solution is then 

(39) h& = £-\a^ - ) + W M p u) e , for some /3« e R K , 

where C~ x is the left pseudoinverse of L M . Here (3^ — (3^(t,x) is a smooth function to be determined. It 
is this arbitrariness in the solution of at order j that allows exactly the freedom necessary to impose the 
solvability condition (38) at order j + 2 of the matching procedure. 

In particular, the leading order (3 of (36) will be determined by the solvability condition (38) at order 2. 
Indeed, at order 1, (37) becomes 

(40) £ M /i (1) =L(vV)M-C (1) (M). 
Its solvability condition (38) is simply 

(41) = (ev • VM) = V- (veM) , 

which is automatically satisfied since the right side vanishes identically when the lattice gas is diffusive (27). 
Solving (40) for M 1 ) then yields 

(42) = C-^Lv • VM - C (1) (M)) + W M f3 (1) e, 
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for some /3^' e H K . At order 2, (37) becomes 

(43) C M hM = [Td t + ±L 2 (vV) 2 ]M + L(vV)hW 

-\V 2 C {0 \M) ■ hWhW - VC^(M) ■ . 

Its solvability condition (38) is 

(44) = (e[ra ( M+iL 2 (vV) 2 M + L(vV)/i (1 ']) 

= Tdt(eM) + LV^veOC" 1 + ±/)v • VM) - iV-(evC (1) (M)) , 

which leads to the evolution equation for [5 = /3(i,x) 

(45) d t (eM) = ^W-(evC (1 \M)) ^V-(((ve(£- 1 + \l)W M ew ■ V) /3) , 

where M = M(f3). The first term on the right is a convection term provided it is nonzero. This will be 
the case whenever M is not a local equilibrium of which can only happen if is not in semidetailed 
balance (recall the H- Theorem). The second term on the right is a diffusion term provided the diffusion 
4-tensor 

(46) ((veTC + \l)W M ev) 

is negative definite. The negativity of ((•ve£~ 1 W M ev) follows directly from Theorem 3.7 and the fact that 
W M ev is not in E since (27) implies that (W M eev) = 0. That this negativity is enough to overcome the 
antidiffusion term in (46) that arises from the second term in the Taylor expansion of the discrete advection 
is a deeper fact due to Henon [12]. This will be made explicit in the specific numerical examples that we 
study later. 

In general, the determination of (3^ at order j + 2 is a consequence of the diffusive property (27) of C 
since, for j > 0, it is seen that a^ +1 ^ has the general form 

(47) a (j+1) = L(v • V)W h J {l) e + , 

where s^ +1 ^ denotes the remaining terms, each of which depends on hf^ through h^ _1 \ but not on (3^\ 
Differentiating the diffusive property (27) with respect to /3 leads to the identity 

(48) (ve W M (/3)/3 w e) = 0® d ? {wSM0)) = . 
The solvability condition (38) at order j + 1 then reduces to 

= (ea (3+1) ) 

(49) = LV-(veW M f3^ e) + (es^ +1 '>) 

= (es^l) . 

This gives a forced, linearization of the convection-diffusion equation (45) that governs the evolution of the 
as yet undetermined /J^' -1 ). In this way one can systematically construct order by order from solutions 
of (45) and its forced linearizations. 
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5 Consistency, Stability and Convergence 

We consider finite truncations of the formal expansion constructed in the last section 
(50) H^(t, X ,p) = J2S k h^(t, X ,p). 

k=0 

By construction, H < - q '(md 2 T, i5L,p) satisfies 

AH^ - - C[Hto] = , 

where formally T^) = 0(S q+1 ). 

Definition 5.1 Let q > be a fixed integer and Bg a finite dimensional Banach space with t\-norm || • 

(i) Consistency 
If 

lim^||^)(m,,.)|U = 0, 

then the lattice Boltzmann method is said to be consistent. 

(ii) Convergence 

If F(m, ■, •) e Bs is the solution to the lattice Boltzmann method (6) and 

A lim o ||F(m,-,-)-ff (9) (m,-,-)||5 = 0, 

for all integers m such that < mAt < T, then the lattice Boltzmann method is said to be convergent. 
(Hi) Stability 

Define the block diagonal matrix L$ : Bs — > Bs where the l-th diagonal block is defined by 



C{F,H^)h= [ \l + VC((1- s)F + sH^) 
Jo L 



• hds , 

for feel. The lattice Boltzmann method is said to be stable if for some r > 0, the family of matrices 

lf[L s :0<At<TandO< mAt < T 1 , 
U=o J 

is uniformly bounded. 

We now prove a theorem that resembles the easy direction of the classical Lax Equivalence Theorem 
found in most finite difference texts, see [21] for example. 

Theorem 5.2 Suppose a lattice Boltzmann method is consistent. Then stability is a sufficient condition for 
convergence. 

Proof: Let F(m, i, •) be a solution to the lattice Boltzmann equation (6) and 

E = F — . 
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Note that 

C(F) - C{H {q) ) + E = C(F, H {q) )E . 

Also, 

C{F) - C{H^) =AE-E + T<«) , 

where 

A lirn o ^||^)(m,,.)||, = 0. 

Hence, 

£(F,H^)E = AE + T( q K 
There exists a permutation matrix P such that 

E(m + 1, •, •) = PLsPtp ■ E{m, •, •) + PT^(m, •, •) . 

Let = PT( q \m,-,-) and Lg = PLgP 1 . Then, by stability and consistency, there exists a constants 
C\ , C2 so that 

m / m— 1 \ 

\\E(m+l,.,.)\\ s = ||E II Ls)T {q Hn,;-) 

n — \r— n+1 / 
m— 1 

< Ci2||T(«)(n,v)||* + ||T(«)(m,.,.)|| 4 

< (m- l)CiC 2 At + C 2 At 

< mCAt . 

Hence, the method is convergent. □ 

We now establish sufficient conditions for stability of a lattice Boltzmann method using the ideas of 
monotone difference methods, e.g., [23]. Consider the operator B defined on I having the p-th coordinate 
function given by 

(51) B[F](p) = F( P )+C[F}(p). 

The derivative of B is the \P\ x \P\ Jacobian matrix Jb(F) whose p, q-th entry is 

We assume Jb{F) to be a continuous function of F. 
Definition 5.3 Let 

fc=0 

be a \P\- dimensional interval upon which Js is a nonnegative matrix. That is, F € £ implies that each entry 
of Jb{F) i s nonnegative. Then £ is called a domain of monotonicity of the lattice Boltzmann method 
(6). The vectors 

M = [M_(p)] p£P and M + = [M + (p)] peP 
are called the extreme points of £ . 
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The following theorem demonstrates the invariance property of the advection operator on a domain of 
monotonicity. 

Theorem 5.4 Let £ be a domain of monotonicity for a lattice Boltzmann method (6) with extreme points 
M and M+. If 



then B leaves £ invariant. That is, F G £ implies BF e £ . 

Proof: Note that Jb > on £ implies the coordinate function B[F] (p) is monotonically increasing. Moreover, 
continuity of Js implies Jg(M + ) = so that M + {p) maximizes B[F](p) on [M_(p), M + (p)]. Hence, 



A stability condition can be established for lattice Boltzmann methods whose collision operators have 
the following conservation property. 

Definition 5.5 A collision operator C is said to conserve mass if (C(F)) = for all F el. That is, e = 1 
is a locally conserved quantity. 

Theorem 5.6 Let £ be any domain of monotonicity for the lattice Boltzmann method (6) such that the 
collision operator is zero at the extreme points and suppose F(0,i, •) and H^(0,i,-) are vectors in £. If C 
conserves mass, then the method is stable. 

Proof: Since Jb{F) > and the fact that \\ ■ \\s is the t\ norm, we get 



where the last equality follows from the fact that C conserves mass. Since by Theorem 5.4 B leaves £ 



C(M_) =C(M+) = 0, 



B[F](p) < B[M + ](p) = M + (p) + C[M+](p) = M+(p) . 



□ 




1 




)-s(H(m,i,-)-F(m,i r ))]e£, 



to, i, •) are both in £ for all integers m. Thus, the connectedness 



0<s< 1. 



Hence for G E X we get 




= \\G\\s. 
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□ 



6 Example - A Nonlinear Diffusive System 

We consider a lattice Boltzmann method, called LB1, constructed from a lattice gas on a periodic square 
lattice. A particle can be in one of p = 0, 1, 2, 3 possible states at a node, as are indicated in Table 1. Here, 



(52) 



v(0) = 



, v(l) = 



v(2) = 



, v(3) = 



The collision rules are illustrated in Table 2 and are formally tabulated in Table 3. Detailed balance is 
verified by examining each collision rule. For example, comparing rule 3 with rule 12 yields 

S((0, 0, 1, 1), (1, 1, 0, 0)) = 5((1, 1, 0, 0), (0, 0, 1, 1)) . 

Scmidetailed balance is verified in the same manner. If F{k) = Fk,k = 0,1,2,3, then the generalized 
Boltzmann collision operator C is given by 

(53) 

C(F)(k) — FkFk+iFk+2Fk+3 + FkFk+iFk+2Fk+3 
+-Ffe-Ffc+i-Ffe+2-Ffe+3 — FkFk+iFk+2Fk+3 
— -Ffe-Ffc+i-Ffe+2-Ffe+3 — FkFk+iFk+2Fk+3 

where the sub-indices are evaluated modulo 4. Note that C conserves mass. 

We know from the 7J-theorem (Theorem 3.6 (i)) that at a local equilibrium, C(F) = . Thus, at a local 
equilibrium we have 

C(F)(0)=C(F)(1) => F = F 1 , 
C(F)(0)=C(F)(2) => F = F 2 , 
C(F)(0)=C(F)(3) => F = F 3 . 

That is, 

(54) F =F 1 =F 2 = F 3 = u. 

If we take partial derivatives of (53) and evaluate them at a local equilibrium (noting that (54) holds), then 
the linearized collision operator is the singular symmetric matrix 



-M 



V 

4 



where v = — 4u(l — u). The eigenvalues of C are 

(Ao, Ai, A 2 , A3) 



= (0, v, v) 
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Table 1: Particle States for LB1 (and LB2). 



Particle State p 





1 


2 


3 



































Table 2: Collision Rules for LB1. Configurations that involve changing particles' directions are marked 



Configuration 


Pre-Collision State 


Post-Collision State 


No particles 


















One particle 


















Two orthogonal 
particles* 


















Two head-on 
particles* 


















Three particles 


















Four particles 
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Table 3: Collision Rules for LB1. 



XVU1C 


n 


ri 


S (n, n! ) 




n(0) 


n(l) 


n(2) 


n(3) 


ra'(0) 


n'(l) 


n'(2) 


n'(3) 































1 


1 











1 











1 


1 


2 








1 











1 





1 


3 








1 


1 


1 


1 








1 


4 





1 











1 








1 


5 





1 





1 


1 





1 





1 


6 





1 


1 





1 








1 


1 


7 





1 

1 


1 


-i 
1 


u 


1 


1 


-1 

1 


1 


8 


1 











1 











1 


9 


1 








1 





1 


1 





1 


10 


1 





1 








1 





1 


1 


11 


1 





1 


1 


l 





1 


1 


1 


12 


1 


1 














1 


1 


1 


13 


1 


1 





1 


l 


1 





1 


1 


14 


1 


1 


1 





l 


1 


1 





1 


15 


1 


1 


1 


1 


l 


1 


1 


1 


1 



and the associated unnormalized eigenmatrix is 



(55) 



Q = [qo , qi , q2 , q3] 



110 1 
10 1-1 

1-10 1 
1 0-1-1 



By Theorem 3.7, E = Null(£) = Span[q ]. The pseudoinverse of C is given by 

(56) C- 1 = v- 1 [§ (q iq * + q 2 q*) + iq 3 q^] . 

Clearly, the lattice gas is diffusive (Definition 4.1). A Hilbert expansion (31) for LB1 is determined where 
we take = C and to be the null operator. It follows from (36) that 

(57) /i<°>=u(t,x)qD. 



From (42), (56), and (57), 
(58) 

for some e R. 

Equation (43) yields 

where 



h W =Lv 1 [u x q 1 +u y q 2 }+ f3 (1) q Q 



Ch& = a (2) - 



>(2) 



[Td t + U 2 (v • V)><°> + L(v ■ V)hW 



r (2) = -^ cih (0)y h (l) h (l) 
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If we substitute (57) and (58) into the solvability condition (44), we get the nonlinear diffusion equation 

(59) d t u = ^-D(u)Vu, £>(„) = _ ^1 + 0, n=^. 

The solvability condition for yields 

(60) d t (3 w = AiV • [d(u)V/3 (1) + D'(u)(3^Vu 

so that we can set = in (58). 
The solvability condition for M 4 ' is 

(61) dtl3 (2) = mV- \d(u)X7i3 ( V + D\u)pW Vul - T . 

where T is a smooth function of u and its derivatives, and the solvability condition for is 

(62) d t p {3) = mV- [l»(m)V/3 (3) +D'(u)/3 (3) Vu . 

We can thus take f3^ = 0. The expressions for ftW, j = 0,1,2,3 are given in Appendix A. 
We now consider the grid function 

3 

H^(m,i,-) = J2h {k \m,i,-)6 k . 



fe=0 



Then, 



where 



AH (3) - # (3) - C(H (3) ) = -T(# (3) ) + higher order terms, 



j=o 

Since T<°) = Ch^ and T« = Ch^ - + j = 1, 2, 3, it follows that T« = 0, j = 0, 1, 2, 3. Moreover, 

T (4) = -E4 4) q. = -4 4) q 3 , 

fc=0 

where 2^ 4) is given in Appendix C. Since u e C 4 ([0, l] 2 , (0, 1)) and (3^ e C 2 (R 2 , R), we have 

and the method is consistent. 

Stability of LB1 will follow from Theorem 5.6 once we have established an invariant domain of mono- 
tonicity for the method. 

Theorem 6.1 The set £ = [M_,M+] 4 , where 

M+ = i(l + l/V5) and M_ = 1(1 - l/VE) , 
is a domain of monotonicity for LB1. Further, the collision operator is zero at the extreme points of E. 
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Proof: The elements of the Jacobian matrix J[B\ are (using (53)) 

J[B]k,k = 1 + F k+2 F k+3 + F k+1 F k+3 + F k+ iF k+2 - F k+1 - F k+2 - F k+3 , 

<J[B]k,k+i = —3F k+2 F k+3 + F k F k+3 + F k F k+2 + F k+2 + F k+3 — F k , 

J[B]k,k+2 = — 3F k+ iF k+3 + F k F k+3 + F k F k+ i + F k+ i + F k+3 — F k , 

J[B]k,k+3 = -3F k+ iF k+2 + F k F k+2 + F k F k+ i + F k+ i + F k+2 — F k , 

where the sub-indices are evaluated modulo 4. Note that J[B]{F) will be nonnegative on £ if each of the 
functions 

f{x,y,z) = 1 + yz + xz + xy — x — y — z , 
g(x,y,z) = -3yz + xz + xy + z + y - x 

are nonnegative for M_ < x,y,z < M + . This will be case if the local and boundary minima of / and g are 
nonnegative. 

If we examine the gradients of these functions, then we see that they are both zero at the point x = y = 
z = i. However, this point is not a local extremum for either / or g since their Hessian matrices D 2 f, D 2 g 
are neither positive nor negative definite there. On the boundaries we have 

0</,5< I- 

Hence, J[B](F) > for F e £. 

Since the extreme points of £ are M + = M + q and M_ = M_q , it follows from (57) that 

C(M+) =C(M_) = 0, 

and the proof is complete. □ 

We can now apply Theorem 5.2 to show that the solution to LB1 converges to the solution u of (59) as 
the lattice spacing is refined. 

Numerical Verification 

We compute solutions to the nonlinear diffusion equation (59) for (x, y) G [0, 1] x [0, 1], £ > 0, /z = ^, periodic 
boundary conditions, and initial condition 

u(x, y, 0) = sin(27rx) sin(27ry) + \ . 

The solutions are computed using both LB1 and an explicit finite difference method that is first order accurate 
in time and second order accurate in space. Note that the initial condition for the lattice Boltzmann method 
lies within the domain of monotonicity of Theorem 6.1. 

The lattice Boltzmann approximations on JV x iV grids (N < 256) at lattice point i and time m arc 
computed by 

{u (N) T= l {Fh 
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Figure 1: LB1: u(t,x,y). 

The finite difference solutions are computed only on a 256 x 256 grid and rendered on coarser grids by 
pointwise projection to yield the solutions (v^)" 1 . Note that the finite difference computed solution retains 
its accuracy when projected onto the coarser grids. 

Figure 1 exhibits the solution for LB1 at times t = and t = 1/32. Table 4 and Figure 2 support the 
theoretical second-order accuracy where denotes the error vector which has components (v^)" 1 — 

( u W^m_ Moreover, the ratio of the errors over grid sizes a factor of two apart yields the expected ratio of 
four. 
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Table 4: LB1: Norm comparisons 



N 


\\E {N % 


t = 1/32 


t = 1/16 


t = 3/32 


t= 1/8 


8 


0.00241 


0.00301 


0.00241 


0.00170 


16 


0.000605 


0.000759 


0.000624 


0.000449 


32 


0.000151 


0.000190 


0.000157 


0.000114 


64 


0.0000366 


0.0000465 


0.0000387 


0.0000279 


128 


0.00000804 


0.0000107 


0.00000902 


0.00000662 




N 


\\EW\\ ti /\\EM\\ ti 


8 


3.991 


3.961 


3.867 


3.776 


16 


4.010 


4.002 


3.987 


3.955 


32 


4.124 


4.078 


4.051 


4.062 


64 


4.550 


4.336 


4.287 


4.223 




Figure 2: LB1: Norm comparisons , (a) Errors; (b) error ratios. 
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Table 5: Collision Rules for LB2. Configurations that involve changing particles' directions are marked 



Configuration 


Pre-Collision State 


Post- Collision State 


No particles 


















One particle* 


















Two orthogonal 
particles* 


— ► 




-* — 












Two head-on 
particles* 


















Three particles 


















Four particles 



















7 Example - Another Nonlinear Diffusive System 

We now consider a variation of LB1, called LB2, where again the lattice gas is defined on a periodic square 
lattice, the possible particle states are in Table 1 and the velocity vectors are given by (52). The collision 
rules for LB 2 are given in Table 5. 

The collision operator for LB2 is in semidetailed balance and is given by 

F k F k+1 F k+2 F k+3 + F k F k+ iF k+2 F k+3 
+F k F k+ iF k+2 F k+3 - F k F k+ iF k+2 F k+3 
—F k F k+ iF k+2 F k+3 — F k F k+1 F k+2 F k+a 
+F k F k+ iF k+2 F k+3 — F k F k+ iF k+2 F k+a , 

where the indices are evaluated modulo 4 and, as before, F k = F(k). Note that the collision operator for 
LB2 is identical to the collision operator for LB1 with the exception of the last two terms. The lattice gas 
is diffusive and the collision operator conserves mass. 
At a local equilibrium of C we have 



C(F){k) 

(63) 



F =F 1 =F 2 =F 3 =u. 
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The linearized collision operator is given by 



■(!-«) 



-(2m -1) u 1 

u — (2u — 1) u 
1 u -(2m +1) 



u 
1 
1 



u 1 u -(2u+l) 

The eigenvalues of £ are 

(Ac, Ai, A 2 , A 3 ) = (0, -2(1 - u 2 ), -2(1 - u 2 ), -4u(l - u)) 

and the unnormalized eigenmatrix is the same as for LB1, cf. (55). Let v = \ 2 = — 2(1 — u 2 ). 

The determination of the Hilbcrt expansion for LB2 proceeds in the same manner as for LB1. In this 
case, we get 

h,W = u(t,x)qo, 
u t = uV • D(u)Vu , D(u) = 
= 0, 



(64) 



1 1 

u + 2 



EL 

2T 



d t (3 {2) = — V- \p(u)V(3 {2) + D'(u)f3 (2) Vu 
/3< 3 > = 0, 



where J 7 is a smooth function of u and its derivatives. Consistency of LB2 follows in the same manner as in 
LB1. Indeed, 

AH (V _ ff(3) _ C(ff ( 3)) = T(H ( 3)) = _g(4) q3 ) 

where c 3 4 ^ is given in Appendix C. Consistency then follows from the fact that T(_ff( 3 )) = 0(<5 4 ). 
Theorem 7.1 TTie sei £ = [M_,M+] 4 7 w/iere 



M+ 



§ and M_ = | , 



is an domain of monotonicity for LB2. Further, the collision operator is zero at the extreme points of £ . 
Proof: The elements of the Jacobian matrix J[B] are (using (63) 

J\JS\k,k = ^fe+2-Ffe+3 + Fk+\Fk+2 — -Ffc+2 , 

J\p\k,k+\ = — 2Ffc + 2-Ffc + 3 + FkFk+2 + -Ffe+3 j 

i7[£>]fc,fc+2 = — 2ffe + iFfc + 3 + FkFk+s + FkFk+i — Fk + 1 , 
i/[£>]fc,fc+3 = — 2ffe+i-Ffe+2 + FkFk+2 + -Ffc+i , 

where the sub-indices are evaluated modulo 4. We show that ,7[Z?] is nonncgativc on £ by showing that each 
of the functions 

f(x,y,z) = yz + xz-y 
9(x,y,z) = -2yz + xy + z 
h(x,y,z) = —2yz + xz + xy — x + l 
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Figure 3: LB2: u(t,x,y). 



are nonnegative on M_ < x,y,z < M + . 

The gradients of these functions are never zero in the domain of restriction for the independent variables 
and so the extrema of these functions is located on the boundaries. A lengthy analysis of these functions 
shows that on these boundaries we have 

0<f,g,h<£, 

cf. [9]. Hence, J[B](F) > for F E £. 

Since the extreme points of £ are M + q and M_Q , it follows C(M+) = C(M_) = and the proof is 
complete. □ 

We can now apply Theorem 5.2 to show that the solution to LB2 converges to the solution u of (64) as 
the lattice spacing is refined. 

Numerical Verification 

We repeat the same experiments as was done for LB1. That is, we solve (64) with the same boundary 
conditions and the initial condition 

u(x, y, 0) = | sin(27rx) sin(27rj/) + | . 
The results are given in Figures 3 and 4 and in Table 6. 
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Table 6: LB2: Norm comparisons for A = 1/12 and B = 3/4. 



N 




E iN) ^ 


t = 1/32 


t = 1/16 


t = 3/32 


t = 1/8 


1 8 


0.000692 


0.00202 


0.00144 


0.000827 


16 


0.000276 


0.000443 


0.000321 


0.000198 


32 


0.0000722 


0.000106 


0.0000782 


0.0000491 


64 


0.0000180 


0.0000264 


0.0000197 


0.0000123 


128 


0.00000463 


0.00000668 


0.00000502 


0.00000316 




N 






8 


2.504 


4.572 


4.480 


4.186 


16 


3.827 


4.166 


4.106 


4.021 


32 


4.009 


4.028 


3.971 


3.995 


64 


3.893 


3.948 


3.924 


3.893 




Figure 4: LB2: Norm comparisons for A = 1/12 and B = 3/4. (a) Errors; (b) error ratios. 
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Particle State p 


1 


-i 


h — 


— H 



Figure 5: Particle States - 1-D Burgers' Equation. 



Prc-Collision 


Post-Collision 


71 


ri 


S(n, n') 






1 










o=i(l + e) 










a=i(l-e) 










a=i(l + e) 






\+ 


1 >- 


a=i(l-e) 










1 







Figure 6: Collision Rules — Burgers' Equation. 

8 Example - Burgers' Equation 

Boghosian and Levermore [3] introduced a lattice Boltzmann method for solving the one-dimensional viscous 
Burgers equation 

(65) p t + pp x = vpxx ■ 

The lattice in this case is one-dimensional and periodic, and the particle states are given in Figure 5. The 
collision rules are listed in Figure 6, where the probability of an advection to the right, i.e., in the direction 
of v = 1, is a = 5(1 + e) and to the left in the direction of v = —1 is a = 1 — a. Here, < e < 1 is given. 
The generalized Boltzmann collision operator is 

(66) C(F) = i[-F 1 +F_ 1 +e(F 1 +F_ 1 -2F 1 F_ 1 )] 

= CW(F) + ^C«(F), C>0. 
Thus, if we assume e = CS, then 

C(F) = C (0) (F) + 5C W (F). 
Clearly, is a Boltzmann collision operator in semidetailed balance. 



+1 

-1 
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We have at an equilibrium 



(67) 
so that 



C (0) OF)(1) =C^ 0) (F){-1) 



Fi = F-i = u. 



The linearized collision operator of is 



2 



-1 +1 
+ 1 -1 



The eigenvalues of C are given by (Ao, Ai) = (0, —1), and the unnormalized eigenmatrix is 

Q = [qo,qi] 

The pseudoinverse of C is 



+ 1 +1 
+ 1 -1 



+1 -1 
-1 +1 



C + = -iqiq! = -\ 
We determine the Hilbert expansion (31) as before. In this case 

= uq , 



and for j = 1,2,3,4 



h U) =/?«q + c 0') qi . 



The are given in Appendix B. Letting 



A(u) = u{l-u), c =^f-' V = > 

we have 

(68) d t u + cd x A{u) = [id xx u , 

= 0, 

d t pW +6cd x (A'(u)f3^) = vd xx pW-T, 

= 0, 

where T is a smooth function of u and its derivatives. 
Consider the grid function 

4 

ff<% a, ■) = £><*>(*,*,■)**. 

k=0 

where the are defined from the Hilbert expansion. Then 

AH& - - C(H^) = T(H^) + higher order terms 
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where 

4 

3=0 

It follows that TW = j = 0, 1, 2, 3. Moreover, T( 4 ) = so that 

AH (4) _ #(4) _ C (#(4)) = Q^Sj 

and the method is consistent. Stability follows from the next theorem. 
Theorem 8.1 The set £ = [0, 1] x [0, 1] is a domain of monotonicity for (66). 

Proof: First note that 



Since < e < 1, we see that J[B] > on £. 
The extreme points of £ are 



\{l + e)-eF 1 
\{l + e)-eF_ 1 
i(l-e) + eF 1 
i(l-e) + eF_ 1 . 



M+ = 



M = 



Clearly, C(M+) = C(M_) = and the theorem is proved. □ 



Numerical Verification 

We compute the solutions to the Burgers equation (68) for x £ [0, 1], t > 0, c = 1, v = 2~ 8 , periodic 
boundary conditions, and initial condition 

u(0, x) = \ cos(27rx) + i . 

The finite difference solution v™ of (68) is computed on a grid of size N = 32768 by solving Burgers' equation 
(65) with the conservative monotone difference scheme 



pT +1 = p" 



" M [(pT+if (pT-i) 2 } + 7^=4 [pT+i 2P? + pT-i] 



4Ax l " ri+u J ' (Ax) 2 

and then applying the transformation p = 1 — 2v. The difference method is first order accurate in time and 
second order accurate in space and has the stability criteria At < (Ax) 2 /(2v). Figure 7 exhibits the initial 
condition and the solution at time t = j. 

The lattice Boltzmann solutions ((u^)™ are computed on grids of size N — 128, 160, 192, 256. Here, u™ — 
\{F). Figure 8 exhibits a comparison of the finite difference- and lattice Boltzmann-computed solutions, 
V(t, x) and U(t,x), respectively, at t = -g. This comparison indicates the importance of the underlying 
assumption that e = 0((Ax) 2 ), which weakens as the grid is coarsened. 

Table 7 lists the £i-norm of the error at time t = j. Varying in the table is the grid size for the lattice 
Boltzmann method. Focusing on the ratio ||£ , ( JV ^||^ 1 /||S( 2JV '||^ 1 , the £i-norm results in the table strongly 
support the theoretical 0(5 2 ) convergence of the lattice Boltzmann method. This point is also illustrated in 
Figure 9. The trend towards the value 4 of this quotient is apparent. 
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Figure 8: Burgers' Equation: Finite Difference and Lattice Boltzmann Solutions at t = j. The Finite 
Difference Solution V(t, x) is the case in which N = 32768 . 
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Table 7: Burgers' Comparison: Norm comparisons. 



N 


\\E (N % 


t = 1/32 


t = 1/16 


t = 3/32 


t = 1/8 


256 


4.254xl0~ 3 


7.709xl0~ 3 


9.170xl0~ 3 


1.042xl0" 2 


512 


9.861xl0- 4 


1.775xl0- 3 


2.102xl0" 3 


2.349xl0" 3 


1024 


2.422xl0~ 4 


4.351xl0~ 4 


5.149xl0~ 4 


5.731xl0" 4 


2048 


6.029xl0~ 5 


1.082xl0~ 4 


1.280xl0" 4 


1.424xl0" 4 


4096 


1.504xl0" 5 


2.702xl0~ 5 


3.194xl0~ 4 


3.551xl0" 5 


8192 


3.755xl0~ 6 


6.741xl0~ 6 


7.967xl0~ 5 


8.865xl0- 6 


16384 


9.296xl0" 7 


1.678xl0" 6 


1.977xl0" 5 


2.207xl0~ 6 




N 




256 


4.314 


4.344 


4.362 


4.436 


512 


4.071 


4.079 


4.082 


4.098 


1024 


4.018 


4.020 


4.022 


4.026 


2048 


4.008 


4.006 


4.009 


4.008 


4096 


4.006 


4.008 


4.009 


4.006 


8192 


4.039 


4.017 


4.031 


4.017 




Figure 9: Burgers' Equation: Norm comparisons, (a) Errors; (b) error ratios. 
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9 Conclusion 

We defined a lattice Boltzmann method as an approximation to an ensembled lattice gas method. The 
concept of semidetailed balance for a lattice Boltzmann collision operator was defined and analyzed. This 
property allowed us to prove an £f-theorem which characterized the equilibria of a Boltzmann collision 
operator. An asymptotic Hilbcrt expansion was constructed about an equilibrium solution of a diffusive 
collision operator. Convergence of a lattice Boltzmann method was established by analyzing the behavior 
of a truncated Hilbert expansion as the perturbation parameter approaches zero. Stability, consistency and 
convergence of a lattice Boltzmann method were defined. Stability and consistency were shown to imply 
convergence. Monotone Boltzmann collision operators were also defined and shown to imply stability. 

Three example lattice Boltzmann methods were analyzed and shown to be consistent and stable. These 
properties allowed us to show that the solutions converged; one to the solution of Burgers' equation and the 
others respectively to the solutions of two nonlinear diffusion equations. Numerical results were presented 
that verified the convergence of each of these methods. 
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A Details of LBl and LB2 Examples 

For LBl we have 

fe=i 

where v = — 4u(l — u) and (using the standard notation V = (d x ,d y ) and V = (d x , —d y )) 

i- (i = 0); 

c (°) = 4 ) = 4°) = o. 



2- (j = i); 



3- (j = 2); 



4. (.7=3); 



For LB2 we have 



Lu x , 



= 0. 



(2) _ „(2) 

L 2 [±D'(u)Vu . Vu - V ■ D(m)Vu] . 



„(2) _ l r 2[l n / 
'3 — 2 



,(3) _ l r 3 



,(3) 



L 3 u xxx + LTu xt + [\L 3 d xx + LTd t ] (z/ _1 u x ) 
+ \L 3 d x [(2j/)- 1 L>'(u)Vm • V« - v- 1 V • D(u)V«] 
+ L/?( 2) +ALu- s {Lu v f 

+ \L 3 D'(u)u x [D'{u)Vu ■ Vit - 2V ■ D(u)Vm] ; 

- ii 3 a y [(2^)- 1 L>'(u)Vw • Vu - v- 1 V • D(u)V«] 
+ L/3( 2 )+4i 3 ^- 3 K) 2 



1 L 3 D'(u)u y [D'(u)Vu ■ Vu - 2V • L>(u)Vw] ; 



4 



4 3) = o. 



fe=i 

where Ai = A2 = v = —2(1 — u 2 ), A 3 = — 4w(l — u) and 
1. (J = 0); 

4 0) = 4 0) = 4 0) = . 
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2- (j = i); 



3- (J = 2); 



4- (J = 2); 

4 3) 



4 1} = 
4 1} = 
4 1} = 



Lu x , 
. 



„( 2 ) 



,(2) _ 



1 t2 
2 L 



[{2u - l)v- 2 Vu ■ Vu - V ■ D(u)Vu] 



= g^ 3 "^ + ^Tuxt + [|d xx + LTd t ] 

+ \L 3 d x v~ x d x u + 2v- 3 L 3 u x (u y ) 2 + L[3 X 2) 

- LD'{u)vl3 {2) u x - L 2 u x V • D(u)X7u + 2^ 1 L 2 (2u - l^V • Vu 
+ L 3 d x [^X^V ■ D(u)Vu + i^- 2 X^ 1 {2u- 1)V- Vu] , 

= \L 3 Uyyy + LTUyf + [\ 1^ Oyy + LT ' 6 t ] ( V ~ ^ U X ) 

+ \l?dyv~ x dyu + 2u~ 3 L 3 u 2 x + Lffl 

+ L 2 u y [V • D{u)Vu + 2u~ 1 L 2 (2u - l)u y V ■ Vu] 

- L 3 d v ■ D{u)Vu + v- 2 \^ 1 {2u- 1)V- Vu] . 
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B Details of the Burgers Equation Example 

For = + c^qi in Burgers' equation, 

1. 

= CA(u) - Lu x ; 



2. 



,(3) 



-\L z u xxx + C {\l?d x + Td t ) A'(u) + LpW _ QCA'{u)f3^ 
+C \{pWf - (Lu x - CA{u)f 



4. 



c ( 4 ) = o. 
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C Further Details of the LB1 and LB2 Examples 

For 



in LB1 we have 



and in LB 2 we have 



t< 4 > = -4 4) q 3 



+ 
+ 



48 

\L 2 



l; !l!J!J!J 



) + jL T(u xxt - Uy y t) + Td t c { 3 



(2) 



V • V/3 (2) + V 



2J2) 



\LldJV-d, 



c (3) 



±L 4 (d XXXl -dyyy) ■ fl^Vu) + ^TV ■ U) 

iLV- 3 £»'(n) {{u x f - {uyf) 

4,, -3 



((u x ) 2 + K) 2 )v • d(u)Vu - 2^< 2 ) (K) 2 - K) 2 ) 



- \L 2 D'{u)v(3^ [D'{u)Vu ■ Vu - 2V ■ £>(u)Vu] 



+ 



\LD'(u)v {c^u x — c^Uy^j 



(2) 



4g-^ Uyyyy) H~ 4-^ ^i^xxt Uyyt) H~ TdfC^ 

+ ±L 4 V 2 4 2) + ±L (M 3) - 9,4 3) ) 

+ ^£ 4 (d OT -9 yyy ) • iz-'Vtj + ±L 2 TV • fltti/^Vu) 

+ 2L(2u - lJiz-^c^V - 4 3) ) - i 2 ^ 3 [2L 2 (2w - lJi/^Ajf 1 ((u x ) 4 - {u y ) 4 ) 

- L 2 v\z\{u x f + K) 2 )V • D(u)Vu - ((u x ) 2 - (u y ) 2 )] 

- 2L 2 (2u - l)Aa [2(2u - l>~ 2 Vu ■ Vu - V ■ D(u)Vu] . 



